Kristin Broms, Neptune & Co.
October 18, 2017
Static maps
3D maps
Interactive maps
Animated maps
Areal data exploration
Additional packages
As you can see, “spatial data” has come to mean much more than creating a map! Visualizing spatial data is an active area of software development, with many exciting packages and functions becoming available on a regular basis.
library(rgdal)
streams <- readOGR(dsn = "../data/R_GIS_Layers/",
layer = "Cove_Drainage_WGS84")## OGR data source with driver: ESRI Shapefile
## Source: "../data/R_GIS_Layers/", layer: "Cove_Drainage_WGS84"
## with 49 features
## It has 9 fields
summary(streams)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x -109.35383 -109.12624
## y 36.46573 36.94164
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## OBJECTID ComID GNIS_ID GNIS_Name
## Min. : 0 Min. : 0 00003435: 2 Cove Wash: 2
## 1st Qu.: 32992 1st Qu.:103664065 NA's :47 NA's :47
## Median : 65593 Median :103665843
## Mean : 68116 Mean :101439114
## 3rd Qu.:105309 3rd Qu.:126660550
## Max. :136421 Max. :126669825
##
## ReachCode FType FCode Source
## 14080204001554: 2 Min. :460 Min. : 0 DOQQ: 1
## 14080204001982: 2 1st Qu.:460 1st Qu.:46003 NHDH:48
## 14080105001404: 1 Median :460 Median :46003
## 14080105001406: 1 Mean :462 Mean :45264
## 14080105001455: 1 3rd Qu.:460 3rd Qu.:46003
## (Other) :41 Max. :558 Max. :55800
## NA's : 1
## Name
## Middle 1 : 2
## Cove Wash : 1
## Cove Wash Middle: 1
## Cove Wash North : 1
## Cove Wash South : 1
## (Other) :17
## NA's :26
plot(streams, col="blue")## note that plots using base R may be distorted.library(ggplot2)
streams_gg <- fortify(streams) ## or use broom::tidy()
# str(streams_gg)
ggplot() +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "blue", size = 1.5)## note that without a coordinates correction, the plot may be distorted.library(ggmap)
devtools::install_github("hadley/ggplot2@v2.2.0")
library(ggplot2)
myLocation <- c(lon = -109.2279, lat = 36.545)
myMap <- get_map(location = myLocation, source = "google",
maptype = "satellite", crop = FALSE, zoom = 13)
g <- ggmap(myMap, darken = c(0.25, "white"))
g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "lightblue", size = 1.5)
ggsave("shp_and_ggmap.png", dpi = 72) library(dplyr)
alum_data <- plot_data %>% filter(Analyte == "ALUMINUM")
g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
color = "lightblue", size = 1.5) +
geom_point(data=alum_data,
aes(x = Longitude, y = Latitude, size = Result, fill = Result),
shape = 21, alpha = 0.9) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
scale_size(range = c(2, 8)) +
guides(size = FALSE) +
coord_equal()
ggsave("static_map.png", dpi = 72) head(alum_data)## Analyte Result Units Latitude Longitude
## 1 ALUMINUM 9120.859 mg/kg 36.59625 -109.1736
## 2 ALUMINUM 7429.408 mg/kg 36.59625 -109.1736
## 3 ALUMINUM 1657.721 mg/kg 36.61342 -109.1356
## 4 ALUMINUM 1709.072 mg/kg 36.61342 -109.1356
## 5 ALUMINUM 5867.922 mg/kg 36.56392 -109.2120
## 6 ALUMINUM 46050.939 mg/kg 36.56392 -109.2120
Common ArcGIS functions, such as unions and intersections, can be done in R. The main advantages are that all manipulations are recorded and reproducible, and code can be automated.
library(scatterplot3d)
with(alum_data, {
lollipop_plot <- scatterplot3d(streams_gg$long, streams_gg$lat,
rep(0, length(streams_gg$long)), # x y and z axis
color = "blue", pch = 16, type = "p",
angle = 120,
scale.y = 1.75,
main = "Example Lollipop plot",
xlab = "Longitude",
ylab = "Latitude",
zlab = "Concentration values",
xlim = c(-109.273, -109.15),
ylim = c(36.50, 36.59),
zlim = range(0, 46000),
grid = TRUE,
box = FALSE)
# add the legend
legend("topright", inset = 0.07,
bty = "n",
title = "Type of Results",
c("High","Low"),
fill = c("#1B9E77", "#D95F02"))
# add the lollipop points
lollipop_plot$points3d(Longitude, Latitude, Result,
col = "#282830",
pch = 21,
bg = ifelse(Result > 10000, "#1B9E77", "#D95F02"),
lwd = 1,
type = "h",
cex = (Result / 50000) + 1)
})lollipop plot
(Generated using a modified version of the scatterplot3d function: https://github.com/USEPA/R-User-Group/tree/master/contributedCode/scatterplot3dMap)
Reference: Beaulieu, et al. 2016. Estimates of reservoir methane emissions based on a spatially balanced probabilistic-survey. Limnology and Oceanography, 61: S27-S40.
library(rgl)
library(dplyr)
plot_3d <- with(alum_data,
plot3d(Longitude, Latitude, log(Result), # x y and z axis
col = ifelse(Result > 10000, "#1B9E77", "#D95F02"), size = 5,
type = "p"))
rglwidget(elementId = "plot3drgl") # to show in presentation## plotly recommends the development version of ggplot2, but will also work with the
## latest version of the package, version 2.2.1
# devtools::install_github('hadley/ggplot2')
# library(ggplot2)
library(devtools)
library(dplyr)
library(plotly)
p <- plot_ly(alum_data,
x = ~Longitude, y = ~Latitude, z = ~Result,
marker = list(color = ~Result,
colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
add_markers() pinteractive_map <- g +
geom_point(data = alum_data,
aes(x = Longitude, y = Latitude, size = Result, fill = Result),
color = "#000000", shape = 21, alpha = 0.8) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728",
guide = "legend") +
scale_size(range = c(2, 8)) +
coord_equal() +
ggtitle("ALUMINUM")ggplotly(interactive_map + theme_void(), filename="plotly")ggmap and plotly are not completely compatible. (ggmap uses ggplot2.2.0, plotly wants the latest development version of ggplot2.)
“ggiraph” also makes ggplot figures interactive. It currently has limited uses but may be expanding.
library(animation)
ani.record(reset = TRUE)
for(a in unique(plot_data$Analyte)) {
plot_data_a <- subset(plot_data, Analyte == a)
animated_maps <- g +
geom_path(data = streams_gg,
aes(x = long, y = lat, group = group),
size = 1.5, color = 'lightblue') +
geom_point(data = plot_data_a,
aes(x = Longitude, y=Latitude, size = Result, fill = Result),
shape = 21, alpha = 0.8) +
scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
scale_size(range = c(2, 8)) +
guides(size = FALSE) +
coord_equal() +
ggtitle(a)
print(animated_maps)
ani.record()
}
oopts = ani.options(interval = 1)
ani.replay()
saveHTML(ani.replay(), img.name = "animation_plot",
htmlfile = "animation_5metals.html", ani.width=600, ani.height=400)Alternate package: library(gganimate)
Alternate package: library(anim.plots)
# Show neighbors on map of subbasins
library(sp)
class(lulc.sp1)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(lulc.sp1, zcol = "PCTHERB", col = "black",
main = "Percent Herb")library(spdep)
## determine who shares a border
ctchmt.nb1 <- poly2nb(lulc.sp1, queen = TRUE)
# put neighbors into list
ctchmt.nbwts.list <- nb2listw(ctchmt.nb1, style = "B")
moran.test(lulc.sp1$PCTHERB, ctchmt.nbwts.list)##
## Moran I test under randomisation
##
## data: lulc.sp1$PCTHERB
## weights: ctchmt.nbwts.list
##
## Moran I statistic standard deviate = 2.9826, p-value = 0.001429
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.212182899 -0.016949153 0.005901866
par(mar=c(5, 5, 2, 2))
moran.plot(lulc.sp1$PCTHERB, listw = ctchmt.nbwts.list,
labels = lulc.df1$PCTHERB,
xlab = "PCTHERB", ylab = "Spatially Lagged PCTHERB",
main = "Moran Plot for PCT HERB")library(magrittr) ## for %<>% function
library(dplyr)
library(Rgraphviz) ## to convert data to necessary data types
library(visNetwork) ## to create the interactive network graph
network_interactive <- function(g, nodes, coords){
## first, switch graph construction to data frames:
nNodes <- length(nodes)
nodes_df <- matrix(NA, nrow=nNodes, 4) # save first 4 attributes
for (i in 1:nNodes){
nodes_df[i, ] <- unlist(jsonlite::fromJSON(nodes[[i]]))[1:4]
}
nodes_df <- as.data.frame(nodes_df, stringsAsFactors=FALSE)
## because edges uses the names, need id to equal names of nodes, not numbers
names(nodes_df)[1:3] <- c("old_id", "id", "type")
## and id column needs to go first: (only for igraph package, not visNetwork)
## igraph packge was used to help decide node coordinates
nodes_df <- nodes_df[, c('id', 'old_id', 'type')]
g_edges <- Rgraphviz::buildEdgeList(g)
nEdges <- length(g_edges)
edges_df <- NULL
for (i in 1:nEdges){
tmp <- c(from=g_edges[[i]]@from,
to=g_edges[[i]]@to,
unlist(g_edges[[i]]@attrs))
edges_df %<>% bind_rows(tmp)
}
edges_df$weight <- as.numeric(edges_df$weight)
## Add layer info/ color column
nodes_df$col_layer <-
ifelse(nodes_df$id %in% grep("MainInput", nodes_df$id, value=T), 1, #
ifelse(nodes_df$id %in% grep("OtherInput", nodes_df$id, value=T), 3,
ifelse(nodes_df$id %in% grep("Midvalue", nodes_df$id, value=T), 4,
ifelse(nodes_df$id %in% grep("MidEqn", nodes_df$id, value=T), 5,
6))))
## add coordinates to determine layout of each node:
nodes_df <- full_join(nodes_df, coords)
## specify colors to use for each col_layer
graph_colors <- c("gold", "darkorange", "tomato",
"palegreen", "seagreen", "royalblue")
# frame = borders
frame_colors <- c("gold3", "darkorange3", "tomato4",
"palegreen3", "seagreen4", "royalblue4")
## (visNetwork doesn't like color names with numbers)
frame_colors_rgb <- rgb(t(col2rgb(frame_colors)), maxColorValue = 255)
## Add attributes so that the graph looks good:
nodes_df$shape <- "ellipse"
nodes_df$color.background <- graph_colors[nodes_df$col_layer]
nodes_df$color.border <- frame_colors_rgb[nodes_df$col_layer]
nodes_df$color.highlight.border <- "darkred"
edges_df$arrows <- "to" # draw arrows on the edges.
edges_df$color.highlight <- "black"
network <- visNetwork::visNetwork(nodes_df, edges_df, width="100%", physics = TRUE) %>%
visNetwork::visIgraphLayout() %>%
visNetwork::visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE,
manipulation = FALSE) %>%
visNetwork::visEdges(color = "slategray", smooth = list(roundess = 1)) %>%
visNetwork::visNodes(font = list(size = 22), shape="ellipse")
return(network)
}network_interactive(g = networkDag, nodes = networkNodes, coords = networkCoords)par(mar=rep(0, 4)) ## change margins so plot fills entire space
## default colors:
plot(1:8, pch=16, cex=5, col=1:8)## Defining a new palette
new_palette <- c("darkred", "chartreuse", "turquoise", "purple",
"gray45", "plum", "black", "#F08800")
palette(new_palette)
plot(1:8, pch=16, cex=5, col=1:8)palette("default") ## return palette to default colorssessionInfo()## R version 3.3.3 (2017-03-06)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] visNetwork_2.0.1 Rgraphviz_2.18.0 graph_1.52.0
## [4] BiocGenerics_0.20.0 magrittr_1.5 spdep_0.6-15
## [7] Matrix_1.2-11 bindrcpp_0.2 plotly_4.7.1
## [10] devtools_1.13.3 dplyr_0.7.3 rgl_0.98.1
## [13] ggplot2_2.2.1.9000 rgdal_1.2-11 sp_1.2-5
##
## loaded via a namespace (and not attached):
## [1] httr_1.3.1 maps_3.2.0 tidyr_0.7.1
## [4] jsonlite_1.5 viridisLite_0.2.0 splines_3.3.3
## [7] gtools_3.5.0 shiny_1.0.5 assertthat_0.2.0
## [10] expm_0.999-2 stats4_3.3.3 yaml_2.1.14
## [13] LearnBayes_2.15 backports_1.1.0 lattice_0.20-35
## [16] glue_1.1.1 digest_0.6.12 colorspace_1.3-2
## [19] htmltools_0.3.6 httpuv_1.3.5 plyr_1.8.4
## [22] pkgconfig_2.0.1 gmodels_2.16.2 purrr_0.2.3
## [25] xtable_1.8-2 scales_0.5.0.9000 gdata_2.18.0
## [28] jpeg_0.1-8 ggmap_2.7 tibble_1.3.4
## [31] withr_2.0.0 lazyeval_0.2.0 proto_1.0.0
## [34] mime_0.5 deldir_0.1-14 memoise_1.1.0
## [37] evaluate_0.10.1 nlme_3.1-131 MASS_7.3-47
## [40] tools_3.3.3 data.table_1.10.4 geosphere_1.5-5
## [43] RgoogleMaps_1.4.1 stringr_1.2.0 munsell_0.4.3
## [46] rlang_0.1.2 rjson_0.2.15 htmlwidgets_0.9
## [49] igraph_1.1.2 crosstalk_1.0.0 bitops_1.0-6
## [52] base64enc_0.1-3 labeling_0.3 rmarkdown_1.6
## [55] boot_1.3-19 gtable_0.2.0 reshape2_1.4.2
## [58] R6_2.2.2 knitr_1.17 bindr_0.1
## [61] rprojroot_1.2 stringi_1.1.5 Rcpp_0.12.12
## [64] mapproj_1.2-5 png_0.1-7 coda_0.19-1